<!DOCTYPE HTML>
<html>
<head>
<meta charset="UTF-8">
<title>Figures 6.21-6.23: Basis pursuit using Gabor functions</title>
<link rel="canonical" href="http://cvxr.com/cvx/examples/cvxbook/Ch06_approx_fitting/html/basispursuit.html">
<link rel="stylesheet" href="../../../examples.css" type="text/css">
</head>
<body>
<div id="header">
<h1>Figures 6.21-6.23: Basis pursuit using Gabor functions</h1>
Jump to:&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#source">Source code</a>&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#output">Text output</a>
&nbsp;&nbsp;&nbsp;&nbsp;
<a href="#plots">Plots</a>
&nbsp;&nbsp;&nbsp;&nbsp;<a href="../../../index.html">Library index</a>
</div>
<div id="content">
<a id="source"></a>
<pre class="codeinput">
<span class="comment">% Section 6.5.4</span>
<span class="comment">% Boyd &amp; Vandenberghe "Convex Optimization"</span>
<span class="comment">% Original by Lieven Vandenberghe</span>
<span class="comment">% Adapted for CVX by Argyris Zymnis - 11/27/2005</span>
<span class="comment">%</span>
<span class="comment">% Here we find a sparse basis for a signal y out of</span>
<span class="comment">% a set of Gabor functions. We do this by solving</span>
<span class="comment">%       minimize  ||A*x-y||_2 + ||x||_1</span>
<span class="comment">%</span>
<span class="comment">% where the columns of A are sampled Gabor functions.</span>
<span class="comment">% We then fix the sparsity pattern obtained and solve</span>
<span class="comment">%       minimize  ||A*x-y||_2</span>
<span class="comment">%</span>
<span class="comment">% NOTE: The file takes a while to run</span>

clear

<span class="comment">% Problem parameters</span>
sigma = 0.05;  <span class="comment">% Size of Gaussian function</span>
Tinv  = 500;   <span class="comment">% Inverse of sample time</span>
Thr   = 0.001; <span class="comment">% Basis signal threshold</span>
kmax  = 30;    <span class="comment">% Number of signals are 2*kmax+1</span>
w0    = 5;     <span class="comment">% Base frequency (w0 * kmax should be 150 for good results)</span>

<span class="comment">% Build sine/cosine basis</span>
fprintf(1,<span class="string">'Building dictionary matrix...'</span>);
<span class="comment">% Gaussian kernels</span>
TK = (Tinv+1)*(2*kmax+1);
t  = (0:Tinv)'/Tinv;
A  = exp(-t.^2/(sigma^2));
ns = nnz(A&gt;=Thr)-1;
A  = A([ns+1:-1:1,2:ns+1],:);
ii = (0:2*ns)';
jj = ones(2*ns+1,1)*(1:Tinv+1);
oT = ones(1,Tinv+1);
A  = sparse(ii(:,oT)+jj,jj,A(:,oT));
A  = A(ns+1:ns+Tinv+1,:);
<span class="comment">% Sine/Cosine basis</span>
k  = [ 0, reshape( [ 1 : kmax ; 1 : kmax ], 1, 2 * kmax ) ];
p  = zeros(1,2*kmax+1); p(3:2:end) = -pi/2;
SC = cos(w0*t*k+ones(Tinv+1,1)*p);
<span class="comment">% Multiply</span>
ii = 1:numel(SC);
jj = rem(ii-1,Tinv+1)+1;
A  = sparse(ii,jj,SC(:)) * A;
A  = reshape(A,Tinv+1,(Tinv+1)*(2*kmax+1));
fprintf(1,<span class="string">'done.\n'</span>);

<span class="comment">% Construct example signal</span>
a = 0.5*sin(t*11)+1;
theta = sin(5*t)*30;
b = a.*sin(theta);

<span class="comment">% Solve the Basis Pursuit problem</span>
disp(<span class="string">'Solving Basis Pursuit problem...'</span>);
tic
cvx_begin
    variable <span class="string">x(30561)</span>
    minimize(sum_square(A*x-b)+norm(x,1))
cvx_end
disp(<span class="string">'done'</span>);
toc

<span class="comment">% Reoptimize problem over nonzero coefficients</span>
p = find(abs(x) &gt; 1e-5);
A2 = A(:,p);
x2 = A2 \ b;

<span class="comment">% Constants</span>
M = 61; <span class="comment">% Number of different Basis signals</span>
sk = 250; <span class="comment">% Index of s = 0.5</span>

<span class="comment">% Plot example basis functions;</span>
<span class="comment">%if (0) % to do this, re-run basispursuit.m to create A</span>
figure(1); clf;
subplot(3,1,1); plot(t,A(:,M*sk+1)); axis([0 1 -1 1]);
title(<span class="string">'Basis function 1'</span>);
subplot(3,1,2); plot(t,A(:,M*sk+31)); axis([0 1 -1 1]);
title(<span class="string">'Basis function 2'</span>);
subplot(3,1,3); plot(t,A(:,M*sk+61)); axis([0 1 -1 1]);
title(<span class="string">'Basis function 3'</span>);
<span class="comment">%print -deps bp-dict_helv.eps</span>

<span class="comment">% Plot reconstructed signal</span>
figure(2); clf;
subplot(2,1,1);
plot(t,A2*x2,<span class="string">'--'</span>,t,b,<span class="string">'-'</span>); axis([0 1 -1.5 1.5]);
xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'y_{hat} and y'</span>);
title(<span class="string">'Original and Reconstructed signals'</span>)
subplot(2,1,2);
plot(t,A2*x2-b); axis([0 1 -0.06 0.06]);
title(<span class="string">'Reconstruction error'</span>)
xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'y - y_{hat}'</span>);
<span class="comment">%print -deps bp-approx_helv.eps</span>

<span class="comment">% Plot frequency plot</span>
figure(3); clf;

subplot(2,1,1);
plot(t,b); xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'y'</span>); axis([0 1 -1.5 1.5]);
title(<span class="string">'Original Signal'</span>)
subplot(2,1,2);
plot(t,150*abs(cos(w0*t)),<span class="string">'--'</span>);
hold <span class="string">on</span>;
<span class="keyword">for</span> k = 1:length(t);
  <span class="keyword">if</span>(abs(x((k-1)*M+1)) &gt; 1e-5), plot(t(k),0,<span class="string">'o'</span>); <span class="keyword">end</span>;
  <span class="keyword">for</span> j = 2:2:kmax*2
    <span class="keyword">if</span>((abs(x((k-1)*M+j)) &gt; 1e-5) | (abs(x((k-1)*M+j+1)) &gt; 1e-5)),
      plot(t(k),w0*j/2,<span class="string">'o'</span>);
    <span class="keyword">end</span>;
  <span class="keyword">end</span>;
<span class="keyword">end</span>;
xlabel(<span class="string">'t'</span>); ylabel(<span class="string">'w'</span>);
title(<span class="string">'Instantaneous frequency'</span>)
hold <span class="string">off</span>;
</pre>
<a id="output"></a>
<pre class="codeoutput">
Building dictionary matrix...done.
Solving Basis Pursuit problem...
 
Calling SDPT3: 61625 variables, 502 equality constraints
------------------------------------------------------------

 num. of constraints = 502
 dim. of socp   var  = 61625,   num. of socp blk  = 30562
*******************************************************************
   SDPT3: Infeasible path-following algorithms
*******************************************************************
 version  predcorr  gam  expon  scale_data
    NT      1      0.000   1        0    
it pstep dstep pinfeas dinfeas  gap      prim-obj      dual-obj    cputime
-------------------------------------------------------------------
 0|0.000|0.000|1.5e+00|1.7e+02|7.6e+06| 4.324221e+04  0.000000e+00| 0:0:02| chol  1  1 
 1|1.000|0.995|3.6e-07|1.0e+00|8.6e+04| 4.270270e+04 -7.964108e+01| 0:0:04| chol  1  1 
 2|1.000|0.989|4.6e-07|2.1e-02|1.1e+04| 1.102337e+04 -7.582646e-01| 0:0:07| chol  1  1 
 3|0.970|0.989|2.1e-08|1.2e-03|3.4e+02| 3.397756e+02  3.356358e-01| 0:0:10| chol  1  1 
 4|0.484|0.886|1.2e-08|2.3e-04|2.1e+02| 2.206064e+02  6.421842e+00| 0:0:14| chol  1  1 
 5|0.772|0.915|2.8e-09|2.8e-05|1.0e+02| 1.106558e+02  1.003458e+01| 0:0:17| chol  1  1 
 6|0.507|0.828|1.4e-09|5.7e-06|6.6e+01| 7.735771e+01  1.171877e+01| 0:0:20| chol  1  1 
 7|0.514|0.678|6.7e-10|1.9e-06|4.0e+01| 5.172854e+01  1.210713e+01| 0:0:23| chol  1  1 
 8|0.362|0.850|4.3e-10|3.0e-07|3.0e+01| 4.235935e+01  1.241274e+01| 0:0:26| chol  1  1 
 9|0.447|0.604|2.4e-10|1.2e-07|2.0e+01| 3.222978e+01  1.253063e+01| 0:0:30| chol  1  1 
10|0.527|0.932|1.1e-10|8.1e-09|1.2e+01| 2.482939e+01  1.268651e+01| 0:0:33| chol  1  1 
11|0.624|0.687|4.2e-11|2.6e-09|6.6e+00| 1.932119e+01  1.274802e+01| 0:0:36| chol  1  1 
12|0.599|0.980|1.7e-11|6.1e-11|3.7e+00| 1.649674e+01  1.281034e+01| 0:0:39| chol  1  1 
13|0.543|0.647|7.7e-12|2.5e-11|2.2e+00| 1.507030e+01  1.282581e+01| 0:0:42| chol  1  1 
14|0.541|0.964|3.5e-12|2.5e-12|1.3e+00| 1.413866e+01  1.283901e+01| 0:0:45| chol  1  1 
15|0.728|0.656|1.1e-12|1.8e-12|5.3e-01| 1.337669e+01  1.284169e+01| 0:0:49| chol  1  1 
16|0.458|0.608|6.3e-13|1.7e-12|3.4e-01| 1.318123e+01  1.284339e+01| 0:0:52| chol  1  1 
17|0.590|0.955|5.3e-13|1.1e-12|1.9e-01| 1.303175e+01  1.284468e+01| 0:0:55| chol  1  1 
18|0.961|0.987|2.5e-12|1.0e-12|4.6e-02| 1.289055e+01  1.284500e+01| 0:0:58| chol  1  1 
19|0.761|0.883|4.2e-12|1.1e-12|1.6e-02| 1.286137e+01  1.284512e+01| 0:1:01| chol  1  1 
20|0.882|0.884|9.5e-12|1.1e-12|2.8e-03| 1.284795e+01  1.284515e+01| 0:1:05| chol  1  1 
21|0.943|0.968|4.8e-11|1.5e-12|2.8e-04| 1.284543e+01  1.284516e+01| 0:1:08| chol  2  2 
22|0.628|0.970|3.8e-11|2.3e-12|1.5e-04| 1.284530e+01  1.284516e+01| 0:1:11| chol  1  2 
23|0.631|1.000|2.2e-11|3.4e-12|7.8e-05| 1.284523e+01  1.284516e+01| 0:1:14| chol  1  2 
24|0.623|1.000|1.1e-11|4.3e-12|4.2e-05| 1.284520e+01  1.284516e+01| 0:1:17| chol  1  2 
25|0.619|1.000|6.8e-12|2.3e-12|2.3e-05| 1.284518e+01  1.284516e+01| 0:1:21| chol  2  2 
26|0.602|1.000|5.0e-12|1.4e-12|1.3e-05| 1.284517e+01  1.284516e+01| 0:1:24| chol  2  2 
27|0.591|1.000|3.4e-12|1.0e-12|7.0e-06| 1.284516e+01  1.284516e+01| 0:1:27| chol  2  1 
28|0.599|1.000|1.1e-11|1.0e-12|3.9e-06| 1.284516e+01  1.284516e+01| 0:1:30| chol  2  2 
29|0.606|1.000|5.3e-12|1.5e-12|2.2e-06| 1.284516e+01  1.284516e+01| 0:1:33| chol  2  2 
30|0.612|1.000|3.0e-12|1.1e-12|1.2e-06| 1.284516e+01  1.284516e+01| 0:1:37| chol  2  2 
31|0.618|1.000|2.4e-12|1.0e-12|6.4e-07| 1.284516e+01  1.284516e+01| 0:1:40| chol  2  2 
32|0.622|1.000|2.1e-12|1.0e-12|3.4e-07| 1.284516e+01  1.284516e+01| 0:1:43|
  stop: max(relative gap, infeasibilities) &lt; 1.49e-08
-------------------------------------------------------------------
 number of iterations   = 32
 primal objective value =  1.28451560e+01
 dual   objective value =  1.28451557e+01
 gap := trace(XZ)       = 3.44e-07
 relative gap           = 1.29e-08
 actual relative gap    = 1.29e-08
 rel. primal infeas     = 2.06e-12
 rel. dual   infeas     = 1.00e-12
 norm(X), norm(y), norm(Z) = 3.5e+00, 9.9e-01, 1.9e+02
 norm(A), norm(b), norm(C) = 6.9e+02, 1.9e+01, 1.8e+02
 Total CPU time (secs)  = 102.93  
 CPU time per iteration = 3.22  
 termination code       =  0
 DIMACS: 1.6e-11  0.0e+00  8.8e-11  0.0e+00  1.3e-08  1.3e-08
-------------------------------------------------------------------
------------------------------------------------------------
Status: Solved
Optimal value (cvx_optval): +12.8452
 
done
Elapsed time is 106.049858 seconds.
</pre>
<a id="plots"></a>
<div id="plotoutput">
<img src="basispursuit__01.png" alt=""> <img src="basispursuit__02.png" alt=""> <img src="basispursuit__03.png" alt=""> 
</div>
</div>
</body>
</html>